PROGRAM NINE_QUAD_TEST
IMPLICIT NONE

REAL(8) :: x(25),y(25)
integer :: ID(6,25),IEN(9,5)
REAL(8) :: a,b,c,d    ! u=0+ax+by, v=0+cx+dy
REAL(8) :: E,v
REAL(8) :: sigx,sigy,txy
REAL(8) :: Fx(8),Fy(8)
integer :: I,J

a=1.0D-5
b=0.D0
c=0.D0
d=0.D0

data x/0.0D0, &
       1.0D0, &
       1.0D0, &
       0.0D0, &
       0.3D0, &
       0.6D0, &
       0.8D0, &
       0.25D0, &
       0.5D0, &   !node 9
       1.0D0, &
       0.5D0, &
       0.0D0, &
       0.2D0, &   !node 13
       0.8D0, &
       0.9D0, &
       0.1D0, &
       0.5D0, &   !node 17
       0.7D0, &
       0.6D0, &
       0.35D0, &
       0.5D0, &   !node 21
       0.85D0, &
       0.6D0, &
       0.2D0, &
       0.5D0/
data y/0.0D0, &
       0.0D0, &
       1.0D0, &
       1.0D0, &
       0.3D0, &
       0.2D0, &
       0.75D0, &
       0.9D0, &
       0.0D0, &     !node 9
       0.5D0, &
       1.0D0, &
       0.5D0, &
       0.15D0, &     !node 13
       0.15D0, &
       0.85D0, &
       0.95D0, &
       0.2D0, &     !node 17
       0.6D0, &
       0.8D0, &
       0.6D0, &
       0.1D0, &     !node 21
       0.5D0, &
       0.9D0, &
       0.4D0, &
       0.5D0/
data((ID(I,J),I=1,6),J=1,25)/1, 1, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             1, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             1, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1, &
                             0, 0, 1, 1, 1, 1/
data ((IEN(I,J),I=1,9),J=1,5)/1, 2, 6, 5, 9, 14, 17, 13, 21, &
                              2, 3, 7, 6, 10, 15, 18, 14, 22, &
                              3, 4, 8, 7, 11, 16, 19, 15, 23, &
                              4, 1, 5, 8, 12, 13, 20, 16, 24, &
                              5, 6, 7, 8, 17, 18, 19, 20, 25/

E=3.0D7
v=0.3D0

sigx=E/(1-v*v)*(a+v*d)
sigy=E/(1-v*v)*(d+v*a)
txy=E/(1+v)/2.D0*(b+c)

Fx(1)=-(txy+sigx)/6.D0
Fy(1)=-(txy+sigy)/6.D0
Fx(2)=-(txy-sigx)/6.D0
Fy(2)=(txy-sigy)/6.D0
Fx(3)=(txy+sigx)/6.D0
Fy(3)=(txy+sigy)/6.D0
Fx(4)=(txy-sigx)/6.D0
Fy(4)=-(txy-sigy)/6.D0
Fx(5)=-txy/3.D0*2.D0           !node 9
Fy(5)=-sigy/3.D0*2.D0
Fx(6)=sigx/3.D0*2.D0           !node 10
Fy(6)=txy/3.D0*2.D0
Fx(7)=txy/3.D0*2.D0            !node 11
Fy(7)=sigy/3.D0*2.D0
Fx(8)=-sigx/3.D0*2.D0          !node 12
Fy(8)=-txy/3.D0*2.D0

OPEN(unit=10 , FILE = "Patch_Test.in", &
               FORM = "FORMATTED", &
               STATUS = "REPLACE", &
               ACCESS = "SEQUENTIAL")

write(10,"('Patch test to 9Q element')")
write(10,"(4(I10))")25,1,1,1
do i=1,25
    write(10,"(7(I10),3(F20.12))")i,(ID(J,I),J=1,6),x(I),y(I),0.0
end do

write(10,"(2(I10))")1,12
write(10,"(2(I10),F20.12,/,2(I10),F20.12)")(I,1,Fx(I),I=2,3)
write(10,"(2(I10),F20.12,/,2(I10),F20.12,/,2(I10),F20.12)")(I+4,1,Fx(I),I=5,7)
write(10,"(2(I10),F20.12,/,2(I10),F20.12)")(I,2,Fy(I),I=2,3)
write(10,"(2(I10),F20.12,/,2(I10),F20.12,/,2(I10),F20.12)")(I+4,2,Fy(I),I=5,7)
write(10,"(2(I10),F20.12)")4,2,Fy(4)
write(10,"(2(I10),F20.12)")12,2,Fy(8)

write(10,"(4(I10))")9,5,1,0
write(10,"(I10,E20.12,F20.12)")1,E,v
do i=1,5
    write(10,"(11(I10))")i,(IEN(J,I),J=1,9),1
end do

CLOSE(10)

END PROGRAM NINE_QUAD_TEST